function [A, L, X, rho, b] = test_matrix(N)

A = full(delsq(numgrid('N', N))); % five point laplacian

% Don't use N too large!
[vt,et] = eig(A);               % Get the first three eigenpairs.
                                % L_2 / L_1 = rho, L_3/L_1 = "b"
et      = diag(et);             % put eigval's in a vector
[e,idx] = sort(et,1,'descend'); % sort them, high to low, and save index 
v       = vt; 
for i = 1:length(e)
    v(:,i) = vt(:, idx(i));     % sort the vectors
end

% Reference values.
L   = e(1);
X   = v(:, 1);
rho = e(2) / e(1);
b   = e(end) / e(1);

end